function [Tthreshold,NumberRejected,j_star] = BH(pvalue,delta) 


M = length(pvalue);
% sort tests from most to least significant
sorted_pvalue = sort(pvalue);
% compute the threshold for each sorted hypothesis
threshold = delta .*((1:M)/M);
% find significant hypotheses
% j star represents the position of the first hypothesis that is rejected in the set of ordered hypotheses
% starting from the lest significant and working up
j_star = NaN;
for j = M:-1:1
    if sorted_pvalue(j) <= threshold(j) 
        j_star = j;
        break
    end
end
if ~isempty(j_star)&&j_star>0
    NumberRejected = j_star;
    % convert pvalue into tstat
    Tthreshold = -norminv(sorted_pvalue(j_star)/2,0,1);
else
    NumberRejected = M;
    % convert pvalue into tstat
    Tthreshold = -norminv(max(sorted_pvalue)/2,0,1);
end
Tthreshold = max(1.96,Tthreshold);